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Abstract. The coefficients of interatomic potential of simple form Exp-6 for 
neon are obtained. Repulsive part is calculated ab-initio in the Hartree-Fock 
approximation using the basis of atomic orbitals orthogonalized exactly on 
different lattice sites. Attractive part is determined empirically using single 
fitting parameter. The potential obtained describes well the equation of state 
and elastic moduli of neon crystal in wide range of interatomic distances and 
it is appropriate for molecular dynamic simulations of high temperature prop- 
erties and phenomena in crystals and liquids. 



1. Introduction 

Investigation of strongly anharmonic nonlinear atomic systems by molecular dy- 
namics (MD) method at high temperatures, pressures, or study of systems effected 
by large amplitude excitations requires high accuracy of interatomic potential (IP). 
Series expansion of the IP in the displacements of atoms from equilibrium posi- 
tions is widely used both in phonon theory and in MD simulation jy], |[. Usually, 
fourth-order anharmonisms or lower-order ones can be taken into account because 
of complexity of expansion coefficients calculation. As an alternative, realistic po- 
tential method is used j[7[ ^2], [H], in which exact equations of motion of atoms 
are solved using IP of concrete substance without series expansion. Owing to that, 
all-order anharmonisms are taken into account automatically. This advantage of 
realistic potential method is especially useful in the MD simulation of soliton solu- 
tions where atoms approach each other closely. Realistic IP should have simplest 
form to reduce calculation expenses as well as it must describe precisely the proper- 
ties of the substance under extreme conditions. The aim of this paper is to obtain 
such IP. 

Conventional way of realistic IP determination is empirical fitting to the prop- 
erties of gas or a crystal near the equilibrium point |?], [|. However, such poten- 
tials become unreliable at small interatomic distances like to that arising in soliton 
waves. The properties of highly compressed matter (e.g., for neon up to 1 Mbar 
|l3| ) could give an information for obtaining all-distance reliable IP. However, the 
set of properties, which can be measured accurately at megabar pressures is re- 
stricted strongly. Practically, the equation of state and bulk modulus only may 

For Cik modules the precision worsens drastically 



be included in this set 1 13 
even at kilobar pressures 19| , |23| . There is insufficiency of empirical information 
for fitting all the parameters of IP, and ab-initio calculation is required. 

Realistic IP via interatomic distance is obtained in present work for the crys- 
tal and dimer of neon. Repulsive part of the potential is calculated ab-initio in 
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Hartree-Fock approximation using the basis of localized atomic orbitals orthogo- 
nalized exactly on different lattice sites. Attractive part is chosen to have standard 
Van-der-Vaals form of Cr~ 6 with single empirical parameter C. Used approxima- 
tions and calculation details are described in the section 2. In the section 3 repulsive 
part of IP is interpolated by exponential function of interatomic distance (Exp-6 
potential) and the IP parameter are determined. Experimental verification of the 
IP obtained is performed in the section 4 using the data concerning equation of 
state [|l3| D and elastic moduli [g, g, g, [u[ 1§ @ of compressed neon. The IP 
calculated is found to be in a good consistence with the experiment in whole range 
of pressure. 



2. Ab-initio calculation of repulsion potential 

In MD simulations by realistic potential method the problem is divided into two 
stages. The former is quantum-mechanical calculation of the IP at electron level, 
with interatomic distance considering as a parameter. The latter is solving equa- 
tions of motion of atoms using the IP obtained. This division is correspondent to 
adiabatic approximation when atoms and electrons motion is described separately 
I 

Since pair collisions of atoms have maximal probability, we concentrate the at- 
tention on the dimer of neon, and define the IP as a cohesive energy of the dimer. 
Three-atom forces can be taken into account as a correction to the two-atom ones 



using incremental expansion |20[. The estimation of 20 shows three-atom force 
contribution to be small. 

In Hartree-Fock approximation short-range repulsive part of IP is expressed 
through one-electron density matrix. We don't use hard core approximation. Re- 
arrangement of all electron shells is allowed as interatomic distance is altered. 
Localized basis of atomic orbitalls orthogonalized exactly (by Lovdin procedure 



15 1) on different lattice sites is used. In this basis one-electron density matrix has 
the form Q 



p(r'|r; {1}) = 2]>> s (r' - l)tf(r - 1) - £ *v(r' - l')^'i¥>:(r - 1)}, 

Is IV 

(2.1) P = I-(I + S)-\ 

where (p s (r— 1) = Is > is wave function of electron in isolated atom (atomic orbital), 
1 and 1' are radius- vectors of lattice sites, s numerates occupied states of the atom, 
P is orthogonalizing matrix, I is unit matrix, S is overlap integral matrix with the 
elements 



Sj£ =< lV|Ls >;1^1', 

(2.2) S^ a = 0;l = l'. 

We expand repulsive part of IP in the terms of small parameter such as the 
largest overlap integral S. Usually, S << 1 in uncompressed crystal, and over- 
lap integrals grow exponentially as interatomic distance is decreased. The IP is 
expressed through the products of elements of orthogonalizing matrix P and two- 
center Slater-Koster integrals. These integrals are atomic obital matrix elements of 
crystal hamiltoinial operators. The order in S for two-center integrals is estimated 
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using the theorem about average value. The elements of matrix P = I — (I + S) 
are expanded in powers of overlap integrals matrix S 



p5=5»', + 0(S 2 ), 
(2.3) P», =-(S 2 )»,+0(S 3 ). 

The elements of matrix P contain high-order terms along with the main ones 
proportional to S and S 2 . 

Using the estimations described above, we expand the repulsive part of IP in 
powers of S 



(2.4) V sr = £ (0) +W 2 +W 4 + W 6 . 

Here £ (0) is the energy of intcratomis interaction if orthogonalizing of neighbor 
atoms orbitals is neglected, W2,W4,Wq are orthogonalizing corrections. Series 
expansion in S begins for them from the second, the third, and the sixth powers 
respectively. Due to the presence of matrix P, orthogonalizing corrections contain 
high-order terms in S along with the main ones. 



In the equation 2.4 



(2.5) 15< > = ]T E ( ls l^n + V a m + V£\U) + U n 

Is m,m^l 



The first term in equation |2.5| consists of two-center integrals. They are atomic or- 
bital matrix elements of electron-ion interaction potential V™, of neutral isolated 
atom potential V™ , of electron-electron exchange interaction potential V e ™ respec- 
tively. The second term is the energy of nucleus- nucleus interaction. Electron-ion 
interaction potential has the form 

(2.6) V™ = V en (r-m) = -Ze 2 /\r-m\. 
Neutral isolated atom potential is 

(2.7) U a m = V a (r m) = V en (r - m) + 2^ < mt\v c \mt >, 

t 

where 

< mt\v c \m.t >= J ip* (r' — ra)v c {r — r')ip t (r' — m)dr' , 

v c (r-r') = e 2 /\r-r'\. 

Action of electron-electron exchange interaction potential on wave function is de- 
fined as 

(2.8) < ls|V™l ls >= - E < ls ' mi l^|ls, mi > . 



In the equation 2.4 orthogonalizing corrections, W2, W4, We, are of the form 
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W 2 = -2^ E P%><*s>\V i : +Vl\\s>- 
(2.9) J2 P^Pw<y,™t'\v c \ms,H>; 

Imss' tt f ,l^m 

W A = E P "' P "'{ 2 < ls',K> c |lt,ls > - < Ls',li> c |Ls,lt >} - 

-2 J] < E ^ + > + 

\ss' m^l 

+2 E { pl ss> P ™ m < ls ', mi> c |mi, Is > +P™}P}? < Is', mt> c |li, ms >} + 
+2 E -P" 1 ^ 1 ! 2 < Is', U'\v c \mt, ms > - < Is', li> c |ms, mt >} + 

hoass'tt' 

(2.10) +4 E P ™ lp /t< i 2 < 1*', tt>c|li, m* > - < Is', li> c |ms, It >}; 
(2.11) W 6 = - E pl ss'Pw m < Is', nrf> c |I«, m< > . 

lmss'ff .l^m 

Since the orthogonalizing corrections grow exponentially as the interatomic dis- 
tance is decreased it is impossible to say what correction may be neglected. It 
should be checked for each substance under consideration. 

Using the method described, we calculate repulsive part of IP, V sr (equation 
2.4), for neon dimer as a function of interatomic distance d. Atomic orbitals from 
Clementi-Roetti set |l(J are used as a basis. Hartree system of atomic units h = 
e = m e = 1 is applied. The calculation shows the terms and W2 in equation 
2.4 to have the same order of magnitude and opposite signs. These terms are found 
to give major contributions to the IP. The W4 correction consists of 0,02 per cent 
of the IP at equilibrium interatomic distance do. Further, the W4 does not exceed 
of 1 per cent of the IP up to d ~ 0.75<io- Finally, at small d, like to that arising in 
soliton waves (for d above 0,6-0.75 do), the W4 becomes about 2-4 per cent of the 
IP. The contribution of Wq to the IP is small negligibly (0.002 per cent) in whole 
range of d under consideration. 

3. Determination of interatomic potential parameters 

We interpolate calculated points V sr (d) by exponential function of interatonic 
distance using least square method by the formula 



V sr (d) = ^ exp(-a(a; - 1)); 
(3.1) x = d/zo 

with two unknown parameters Aq and a. Experimental equilibrium interatomic 
distance for neon dimer zq — 5,8411 a.u. p8| is used as the third parameter 
of the IP. The parameters are found to be A = (1, 1384 ± 0,0002) • 10~ 4 a.u., 
a = 13, 6407 ± 0, 0037. Interpolation error is 4-1 per cent of V sr when the d is 
altered from equilibrium one to 0.6z . 
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Adding the attractive part, we express the IP in standard Exp-6 form 

V(d) = A exp(-a(a; - 1)) - C(T 6 ; 
(3.2) x = d/z . 

A single unknown parameter C remains in attractive part of IP. We propose to fit 
the C to experimental equilibrium interatomic distance. Using of equilibrium data is 
considered to be reliable at all interatomic distances since the attraction is essential 
near the equilibrium only while ab-initio calculated repulsive part dominates at 
small d. 

For MD simulation of lattice dynamics, it is possible to fit the C to experimen- 
tal data for dimer at T=0 K because the temperature effects will be taken into 
account explicitly, at the stage of equations of motion solving. In this case, for 
neon C=10,7293 (experimental equilibrium interatomis distance in the dimer is 
zo=5,8411 a.u. JI3]). Calculated cohesive energy of dimer is E co h = —1, 4497 • 10~ 4 
a.u., experimental one is E co h = — 1, 338 • 10~ 4 a.u. The discrepancy is 7 per 

cent of experimental value. 

For calculating static properties of a crystal at finite temperature, e.g., equation 
of state, elastic modules, it is better to fit the C to experimental data for a crystal 
at the same temperature. Such determination allows one to take into account 
implicitly three-atom forces, temperature effects, zero-point oscillations, and other 
effects omitted at the stage of IP calculating. In this case, for neon C=7,4030 
(experimental equilibrium interatomis distance in the crystal is do=5,9647 a.u. at 
T = 4,25 K B) . Calculated cohesive energy of uncompressed crystal is E co h — 
—6, 7620 • 10 a.u. per atom, experimental one is E co h = —(7, 35 ± 0, 03) • 10~ 4 
a.u. p6[ . The discrepancy is 7.6 per cent of experimental value. 

4. Results and discussion 
Interatomic potential of neon is given in the figure [l] as a function of inter- 



atomic distance d. The IP calculated by equation 3.2 for dimer is plotted by solid 
curve. Van-der-Vaals constant (C=10,7293) is fitted to experimental equilibrium 
interatomic distance in dimer fug ]. 

"Experimental" IP obtained in |l3| is denoted by solid circles. This IP had been 
determined by interpolating experimental data p(V) (measured at 300 K) by the 
formula Exp-6. The interpolation had been performed in theoretical model taking 
thermal pressure and zero-point oscillations into account explicitl y, e xcluding them 
from the definition of IP. It allows us to compare the 300 K data of |13|] with our zero- 
temperature result. Three-atom forces didn't include explicitly in the model of 
However, in fisjj , the effect of these forces is taken into account implicitly through 
fitting the IP to experimental data for a crystal. In our calculation three-atom 
forces are omitted because of fitting to dimer data. The agreement of calculated 
IP and experimental one indicates that three-atom forces in neon are small at the 
pressures up to IMbar. 

Two remaining curves in the figure ^ are interatomic potentials of neon obtained 
by fitting to experimental data using Lennard-Jones potential (6-12 formula) 

V(x) = e(-2/x 6 + l/x 12 ) 
x = d/z , 
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where e and zq are fitting parameters. Dashed curve is the IP obtained using 
corresponding-states law fitted to vapor-pressure ratio of isotopic liquid |t]] . Dashed- 
dotted curve is the IP fitted to experimental lattice constant and cohesive energy of 
crystal neon at p = 0, T = OK Q. Fitting to equilibrium crystal properties leads to 
bad describing the IP for compressed crystal. Fitting to compressed gas properties 
gives the values of the IP close to experimental ones at moderated pressures. 

Using the IP obtained (eq. |3.2| ) we calculate the equation of state p{V) for solid 
neon. Calculated pressure p against fractional volume is given in the figure ^| as solid 
curve. Van-der-Vaals constant is fitted to experimental equilibrium interatomic 
distance do=5,9647 a.u. measured for crystal neon at T=4.25 K, p~0 Q. 

Experimental points p(V) from |0| (T=300 K) and g] (T=4.2 K) are also given 
in the figure g. At the pressures below 20 kbar theoretical curve is in a good 
agreement with the experimental points of p|. At moderated pressures theoretical 
curve deviates from experimental points of [fj"3[ by 4 per cent. This deviation 
caused, mainly, by neglecting of thermal pressure in our calculation. Figure 
shows temperature sensitivity of the equation of state to be small. 

We calculate bulk modulus of solid neon by means of the IP obtained. Van- 
der-Vaals constant is fitted to experimental equilibrium interatomic distance in the 
crystal Q. Calculated bulk modulus B via the pressure p is given in the figure 
H as solid curve. Experimental points obtained in M at T=4.2 K are plotted as 
solid symbols. Bulk modulus is seen to be more sensitive to the approximations 
used. Growing when the p is enhanced, the difference between calculated B and 
measured one becomes about 7 per cent of experimental B at p=20 kbar. Incorrect 
taking three-atom forces into account at moderated pressures is seems to contribute 
mainly in this discrepancy. In our calculation three-atom forces (and zero-point 
oscillations too) are taken into account implicitly, by fitting the IP to experimental 
data for uncompressed crystal. Thus, calculated B agrees with experimental one 
at small pressures only (to 8 kbar). One can't determine correctly the dynamics 
of alteration of three-atom forces with enhancing of pressure. It is the cause of 
growing the deviation of calculated B from measured one. 

We calculate elastic modules Cik using the IP obtained with Van-der-Waals 
constant fitted to crystal experimental data H. Calculated modules and exper- 
imental ones are given in table [I] for uncompressed solid neon at low tempera- 
tures. Isothermic modules had been obtained in static measurements |^|, pj. Adia- 
batic modules had been measured in ultrasonic and neutron scattering experiments 
|5| H |TJ, |25|, ^2|. However, the difference between isotermic modules and adiabatic 
ones is negligible at the temperatures under consideration (see, e.g., ||). 

The Cik modules are seen to be more sensitive to the measurement method and 
calculation approximations. The difference between theoretical and experimental 
values of C\\ and C44 is about 10 per cent of experimental values for most accurate 
experiment p5f . The agreement is better for C12 modulus (the discrepancy is 
about 2 per cent |25|). The deviation from Cauchy relation S — (C44 — Cyi)ICy2 is 
also given in table[j]. Cauchy violation is the measure of deviation of the IP from 
spherical symmetry. The S = 0, 11 ± 0,03 in p0[, while it falls into experimental 
error bar in other experiments listed in the table III Cauchy relation takes place for 
our calculation results because spherical symmetry form of the IP is supposed in 
theoretical model. Small value of experimental d indicates that spherical symmetry 
approximation for IP is valid for uncompressed neon at least. For another rare 
gas crystal, krypton, experiment j2^] shows Cauchy relation to satisfy well under 
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Table 1. Elastic modules of solid neon 



Rcf. 


T, K 


B, kbar 


C\\, kbar 


C12, kbar 


C44, kbar 


^ (C44— C12 ) 


Method 




4 

2 




4,25 
4,2 


41,12±0,12 
11,0±0,1 










Static mea- 
surements p(V) 




5 
3 




4 
5 


41,36±0,26 
41,2±0,5 










Ultrasonic ve- 
locity measur. 


I 
I 
I 


14 
21 
12 




4,7 
5 
6 


12,1±0,4 
41,24±0,17 
14,52±0,3 


46,9±0,5 
46,61±0,17 
16,49±0,3 


9,7±0,4 
8,55±0,24 
9,03±0,3 


40,0±0,3 
9,52±0,05 
9,28±0,08 


0.03±0,07 
0.11±0,03 
0.03±0,04 


Inelastic neu- 
tron-phonon 
scattering 


Calc. 





40,76 


44,95 


8,67 


8,67 





Ab-initio calc. 



pressure up to 80 kbar. Moreover, for MgO the Cauchy violation is measured to 
drop with enhancing pressure up to 200 kbar f2~4j| . 

Unlike to C%k modules, bulk modulus B is less sensitive to measurement method 
and calculation approximations. The discrepancy of theoretical result and experi- 
mental one doesn't exceed of 4 per cent and falls into experimental error frames. 

5. Conclusion 

Coefficients of realistic IP of simple form Exp-6 are obtained for neon by ab- 
initio calculation of repulsive part in Hartree-Fock approximation in the basis of 
atomic orbitals orthogonalized exactly on different lattice sites. Attractive part is 
determined empirically using single fitting parameter, Van-der-Vaals constant C. 
For fitting the C it is enough to know experimental equilibrium interatomic distance 
in crystal (or dimer), i.e. high pressure experimental data is not required. The IP 
calculated is suitable for molecular dynamic simulations of high temperature and 
high pressure properties and phenomena in crystals and liquids due to simplicity 
of the form and precise describing experimental data in wide range of interatomic 
distances. 
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4,0 4*5 5,0 5,5 6,0 6,5 7*0 

d, a.u. 



Figure 1 . Calculated IP and three potentials fitted to experimen- 
tal data for neon (from ffi, 13, |) respectively). 
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0,0 0,1 0,2 0,3 0,4 0,5 0,6 0,7 



VA/ 





Figure 2 . Equation of state for solid neon (experiment from Q 
(4.2 K) and |l| (300 K)). 
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Figure 3. Bulk modulus of solid neon (experiment from B). 



